Differential gene expression profile of multinodular goiter

Introduction The goiter, a neglected heterogeneous molecular disease, remains a major indication for thyroidectomies in its endemic regions. Objectives This study analyzed differential gene expression in surgical specimens diagnosed with multi nodular and compared the data to that of thyroid tissue without multinodular goiter from patients undergoing thyroidectomy in Manaus-AM, Brazil using RNA-seq technology. Methodology The transcriptome information of the surgical specimen fragments with and without multinodular goiter was accessed by Illumina HiSeq 2000 New Generation Sequencing (NGS) using the RNA-seq NEBNext® Ultra™ RNA Library Prep Kit for Illumina®—#E7530L protocol and differential gene expression analysis. Results Differences were found between the gene expression profiles of the diseased tissues and those of the healthy control tissues; at least 70 genes were differentially expressed. The HOTS gene was expressed only in multinodular goiter tissues (p < 0.05). Conclusion These results demonstrate that the gene expression profile of multinodular goiter is pro-tumoral and that HOTS can play a central role in multinodular goiter development.

Introduction Partial or total removal of the thyroid gland affected by goiter is one of the most commonly performed surgeries in medical practice. The role of goiter as a risk factor for well-differentiated thyroid carcinoma is unclear; however, the prevalence of incidental carcinoma in patients operated for goiter in endemic areas is 10-12% [1,2], which is greater than the overall prevalence of the disease (5.1%) [3].
Although goiter is the main indication for thyroidectomy in goitrogenic geographic areas, its molecular-genetic component has been scarcely studied compared with that of thyroid carcinoma. In addition, studies conducting massive sequencing for thyroid nodular diseases have focused on well-differentiated and undifferentiated thyroid carcinoma [4,5].
Transcriptome analysis, which is gaining prevalence in studies on tumor diseases, allows a better understanding of gene expression profiles in tissues under different conditions, including the knowledge of non-coding RNAs (ncRNAs), monoallelic expression of imprinted genes, and several transcriptional phenomena, such as fluctuations in the expression of nonconstitutive sequences [6].
The performance of ncRNAs, such as the products of the H19 gene, and their relationship with several types of cancer are well reported in the literature. This gene, which is never expressed beyond the embryonic period in normal conditions, has high expression in tumors related to tissue hypoxia and cancer. Aberrant expression patterns of this sequence occur in breast cancer [7] and melanoma [8].
In lung neoplasms, high expression of H19 is related to the epithelial-mesenchymal transition [9]. Its action on metabolic and cell cycle pathways is thought to be involved in the modulation of a pro-tumor state [10][11][12].
The need for a preoperative diagnosis due to gaps in the Bethesda cytological classification from fine needle aspiration (FNA) of thyroid nodules and the advent of molecular studies of these diseases have allowed the development of molecular tests for well-differentiated thyroid carcinoma, notably based on the identification of BRAF and RAS mutations as well as RET/ PTC and PAX8/PPARy rearrangements [13,14], among others such as Afirma GEC 1 , Thy-GenX TEST 1 , and ThyroSeq TEST 1 , all without relevant application for multinodular goiter.
This unprecedented study presents the occurrence of differentially expressed genes between tissues affected by multinodular goiter and disease-free tissues (hereafter referred to as controls) from specimens collected in a geographical area (Amazonas, Brazil) endemic for the disease.

Method
This study was approved by the Human Research Ethics Committee of the Adriano Jorge Hospital Foundation under CAAE 16463813.9.0000.0007 on June 1, 2013.
The study included transcriptome sequencing of two thyroid tissue fragments with multinodular goiter and one control tissue fragment from patients operated in a multinodular goiter endemic region (Manaus, AM, Brazil).
The thyroid fragments used in this study each measured 1 cm 3 and were obtained directly from the surgical specimen after its removal from the cervical region by thyroidectomy. Tissue in the control group was obtained from thyroid tissue fragments from patients with thyroid adenoma, from a region of the thyroid gland 1.5 cm away from the nodule. Tissues were confirmed disease-free by pathology service analysis afterwards.
Immediately after collection, the fragments were stored in microcentrifuge tubes containing the preservative RNAlater™ Stabilization Solution (Thermo Fisher) in a -80˚C deep freezer until histological classification of the specimen by a pathologist as disease-free tissue or tissue with multinodular goiter.
The total RNA was prepared with TRIzol 1 Reagent (lifetechnologies™) protocol, following the manufacturer's recommendations. All following steps to transcriptome sequencing was performed by GenOne Soluções em Biointecnologia Facility (Rio de Janeiro, Brazil). RNA libraries were validated in an Agilent 2100 Bioanalyzer using the RNA 6000 nano Assay. The cDNA libraries were constructed by using NEBNext 1 Ultra™ RNA Library Prep Kit for Illumina 1 -#E7530L RNA-seq protocol with an expected output of 20 GB of data per sample and sequenced in the Illumina HiSeq 2000 platform.

Data analysis
The sequences exploratory analysis were carried out by the Bioinformatics group from the Central Laboratory of High Performance Technologies in Life Sciences (LaCTAD), State University of Campinas (UNICAMP, SP, Brazil). SRA data accession number: PRJNA810866.

Results
The differential expression analysis of two tissue fragments with multinodular goiter and a control tissue identified 65 differentially expressed genes and five pseudogenes, of which 61 were down-regulated and nine up-regulated in thyroid tissue with multinodular goiter compared with the control tissue (Fig 1).
A biological interaction analysis between functional genes identified 423 possible interactions using the Gene MANIA tool with co-expression greater than 50% Fig 2.

Discussion
Medical publications on thyroid surgical diseases are focused on the search for thyroid carcinoma biomarkers [23][24][25][26]. Initial immunohistochemistry and microarray studies comparing the expression profiles of normal, multinodular goiter, adenoma, and carcinoma tissue samples identified different patterns between the diseases but similarity between the groups of genes in the tissue with multinodular goiter and that with papillary carcinoma, which would explain the higher prevalence of incidental carcinoma (preoperatively unknown) in thyroids operated for multinodular goiter in goitrogenic areas and the existence of a common initial tumorigenesis factor [27,28].
The literature is not clear about the molecular origin of multinodular goiter, which certainly involves epigenetic factors, heredity, and the classical iodine deficiency as well as iron and selenium deficiencies in the diet and exposure to foods rich in flavonoids and cyanogenic substances, such as cassava. When chronic, these conditions would lead to mutations and the onset of nodules in the gland [29,30].
In medical practice, when facing thyroid nodules, the presence of malignant lesions needs to be considered [31,32] along with the preoperative FNA investigation and Bethesda's cytological classification, which often need to be repeated, present variable sensitivity and agreement with histopathology, and are inconclusive in up to 30% of cases [33,34].
The need to identify which patients with thyroid nodule should undergo surgery, new therapy strategies, or clinical follow-up justifies the investigation of the molecular characteristics of lesions to determine the risk of multinodular goiter malignancy.
In this study, 70 sequences were differentially expressed between multinodular goiter and disease-free tissue. The down-regulated genes in multinodular goiter were related to several molecular pathways, especially phospholipase C (PLCD4), apoptosis pathways (TNFRSF19), heat shock proteins (HSPA1A, HSPA6), growth factors (SHC3, NRG1), p53 proto-oncogene pathways (THBS1), and chaperone cell repair pathways (BAG3); on the other hand, the inflammatory (COL14A) and complement system (C4B) pathways were up-regulated in multinodular goiter tissue, in addition to the exclusive presence of an antisense transcription from the H19 locus, which encodes the nucleolar protein HOTS in multinodular goiter [35].
These findings are similar to the characteristics of tumor diseases with reduced apoptosis and cellular repair systems along with increased inflammatory activity in the presence of protumor locus products, in this case the HOTS nucleolar protein, which, together with the lncRNA H19, would be possible inducers of cancerous breast, thyroid, liver, kidney, and lung lesions [36,37].
The action of the H19 gene and its products in tumor onset and hyperplastic lesions is evident, with the antagonism of H19 lncRNA and p53 and the activity of one of its gene products, miR-675, in promoting cellular and chromosomal instability being well described, as well as its  hyperexpression in the presence of external factors such as hypoxia [37]. A balance between the products of sense (lncRNA) and antisense (HOTS) H19 transcripts may be related to the regulation of cellular homeostasis.
There have been no studies on multinodular goiter using NGS and RNA-Seq with results similar to those described in the present study. In thyroid carcinoma, some lncRNA are discussed in the gene regulation of disease progression, such as PTCSC3 [38], with XLOC 051122 and XLOC 006074 [39] in local metastasis and PANDAR as a possible target in pro-apoptotic therapies for carcinoma [40].
The question of whether the presence of multinodular goiter can be considered a risk factor for thyroid carcinoma still raises discussion. Recent findings have shown that the same histopathologically diagnosed papillary lesion exhibits different protein expression behavior if the patient has a history of multinodular goiter prior to the diagnosis of neoplasia [40][41][42].
Other studies have shown the importance of membrane proteins in the development of hyperplastic and neoplastic thyroid diseases, especially connexins and aquaporins [43][44][45]. This study identified no differences between the expression profiles of connexins or aquaporins in different tissues, but STARD9 (apolipoprotein) and CPNE4 membrane proteins were down-regulated in multinodular goiter.
Thus, it was possible to identify molecular characteristics of multinodular goiter similar to those found in the genesis of neoplastic tumor lesions, including: 1) reduced cell repair activity; 2) reduced apoptotic pathway activity; 3) increased inflammatory activity; and 4) H19 gene expression with possible inhibitory activity of the p53 proto-oncogene. The presence of H19 gene products hyper-expressed in multinodular goiter, a non-malignant disease with different forms of presentation in endemic regions (small and large multinodular goiters), contributes to the understanding of the genesis of multinodular goiter and its possible roles as a risk factor for malignant lesions and as a possible molecular marker.

PLOS ONE
Further studies in endemic areas with more replicates for NGS analysis and a better understanding of the function of ncRNAs in the development of the disease will be necessary to confirm the hypothesis of multinodular goiter as a pro-tumor state of the thyroid.
Previous findings in the literature have described the low expression of the H19 gene and its products in thyroid cancer [46], which contradicts the high expression of this sequence in the multinodular goiter samples used in this study. This suggests that high H19 gene expression may be used in conjunction with other molecular markers as a diagnostic tool in deciding between conservative and/or surgical treatment for multinodular goiter patients in endemic areas, such as the Amazon.
Future studies should further elucidate the molecular profile of multinodular goiter, deepen the understanding of the functions of non-coding RNA in malignant and non-malignant nodular diseases, and facilitate the development of rapid and cost-effective diagnostic protocols that consider the level of H19 gene expression in patients with multinodular goiter.